TK!2
=v
#1
:n

:s
b
:c
    >>>>  DEFINITE INTEGRALS  <<<<
#3
:n

:s
b
:c
Integral of sin(x) from 0 to pi()/2
#4
:n

:s
b
:c
(exact value is 1):
#5
:n
fun
:v
'sin
:s
i
:c
name of function to be integrated
#6
:n
a
:v
0
:s
i
:u
rad
:d
rad
:c
lower integration limit
#7
:n
b
:v
1.570796326794897
:s
i
:u
rad
:d
rad
:c
upper integration limit
#8
:n
n
:v
128
:s
i
:c
# of intervals (must be even)
#9
:n
res1
:v
1.000000000126001
:s
o
:c
Simpson rule with set # of intervals
#10
:n

:s
b
:c
Integration with precision control:
#11
:n
res11
:v
1.000000000126001
:s
o
:c
  Simpson, iterative procedure
#12
:n
res12
:v
1.000000000126001
:s
o
:c
  Simpson, recursive procedure
#13
:n
res13
:v
1.000000000000135
:s
o
:c
  Romberg, iterative procedure
#15
:n
k
:v
.5
:s
i
:c
parameter for rule function "elinfk"
#16
:n

:s
b
:c
  (left over from model SIMPSON)
=u
#1
:f
rad
:m
57.2957795129
:t
deg
=l
=f
#1
:n
Simpson
:c
Definite integral, Simpson rule
:t
Procedure
:i
fun,a,b,n
:o
value
#2
:n
defint1
:c
Definite integral, Simpson rule with precision control
:t
Procedure
:i
fun,x0,x1
:o
z
#3
:n
defint2
:c
Definite integral, main part of recursive procedure
:t
Rule
:a
fun,x0,x1
:v
z
#4
:n
defint2R
:c
Definite integral, recursive part
:t
Procedure
:i
fun,z0,s0,h,x0,n,eps
:o
z
#5
:n
defint3
:c
Definite integral, Romberg's method
:t
Procedure
:i
fun,a,b
:o
z
#6
:n
elinfk
:c
Integrand for elliptical int. of 1st kind
:t
Rule
:p
k
:a
phi
:v
z
#1
:n
Simpson
:z
" Notation:  fun    name of rule, procedure or built-in function 
"                   defining the integrand
"            a,b    lower and upper limit
"            n      number of integration steps (must be even number)
"            value  value of definite integral

" Description: Standard method of numerical integration using second-degree
"              polynomial approximation of integrand.

if mod(n,2)<>0 then error:= odd_number_of_intervals
h:= (b-a)/n
x:= a
k:= 1
value:= 0
for i=2 to n
    x:= x + h
    k:= 3 - k
    value:= value + k*apply(fun,x)
next i
value:= (2*value+apply(fun,a)+apply(fun,b))*h/3
/E
E
#2
:n
defint1
:z
" Notation:  fun  name of a rule, procedure or built-in function defining
"                 the integrand
"            x0   lower limit
"            x1   upper limit

" This procedure uses Simpson's rule for evaluating definite integrals.  It
" keeps doubling the number of intervals until the difference between the
" current and previous result falls below a defined tolerance limit, which is
" presently set to 1e-8 in the conditional statements.  (Note that either an
" absolute or relative tolerance may be chosen by cancelling and uncancelling
" the appropriate conditional statements.)

" The recursive version of this procedure is called "defintR".  A Simpson
" integrator with a settable number of intervals is in "Simpson".

h:= x1 - x0
s:= apply(fun,x0+h/2)
z:= (apply(fun,x0)+apply(fun,x1))/2 + 2*s
n:= 1
L: z0:= z
   z:= z - s
   h:= h/2
   x:= x0 - h/2
   n:= 2*n
   s:= 0
   for i=1 to n
      x:= x + h
      s:= s + apply(fun,x)
   next i
   z:= z + 2*s
   if abs(z-2*z0)*h/3>1e-8 goto L     " absolute tolerance
   if abs(1-2*z0/z)*h/3>1e-8 goto L   " relative tolerance
z:= z*h/3
/E
31
E
#3
:n
defint2
:z
" This is the main part of a recursive version of a procedure that evaluates
" a definite integral with required absolute or relative tolerance (set to 
" 1e-8 in the last rule).  The iterative version is called "defint".

" Notation:  fun  name of a rule, procedure or built-in function defining
"                 the integrand
"            x0   lower limit
"            x1   upper limit

" Suporting function: defint2R

h = x1 - x0
s = apply(fun,x0+h/2)
z0 = (apply(fun,x0)+apply(fun,x1))/2 + 2*s
z = defint2R(fun,z0*h/3,(z0-s)/3,h/2,x0-h/4,2,1e-8)
/E
E
#4
:n
defint2R
:z
s:= 0
x:= x0
for i=1 to n
   x:= x + h
   s:= s + apply(fun,x)
next i
z:= h*(s0+s/1.5)
if abs(z-z0)>eps then z:= defint2R(fun,z,z/h-s/3,h/2,x0+h/4,2*n,eps)

" Note: This function refers to itself, repeatedly refining the value of the
" integral until the difference between the last two results falls below the
" required absolute tolerance "eps".  If relative tolerance is preferred, the
" above conditional statement should be replaced by the following one:

if abs(1-z0/z)>eps then z:= defint2R(fun,z,z/h-s/3,h/2,x0+h/4,2*n,eps)
/E
14
E
#5
:n
defint3
:z
" Notation:  fun  name of a rule, procedure or built-in function defining
"                 the integrand
"            x0   lower limit
"            x1   upper limit

" This procedure uses Romberg's method for evaluating definite integrals.  
" It uses the trapezoidal rule with doubling number of intervals, and it
" applies Neville's polynomial extrapolation algorithm to extrapolate the
" intermediate approximations to zero interval h.  The process continues
" until the difference between the current and previous results falls below
" a defined tolerance limit.  Presently set to 1e-8, an absolute or relative 
" tolerance may be used.

" Romberg's method is more efficient than Simpson's rule, especially when
" evaluation of the integrand requires extensive computation.

" The polynomial values for Neville's scheme are stored in the temporary
" list "aux".  The trapezoidal rule approximations are in the elements
" 1,2,4,...,(i^2-i)/2+1 .  The Romberg approximations are in the elements
" 1,3,6,...,i*(i+1)/2 .

h:= (b-a)
'aux[1]:= h*(apply(fun,a)+apply(fun,b))/2
i:= 1
k:= 1
jj:= 1
z:= 'aux[1]
L: z0:= z
   x:= a - h/2
   Int:= 0
   for j=1 to jj 
      x:= x + h
      Int:= Int + apply(fun,x)
   next j
   h:= h/2
   jj:= 2*jj
   k:= k + 1
   'aux[k]:= 'aux[k-i]/2 + Int*h
   i:= i + 1
   c:= 1
   for j=2 to i
   c:= 2*c
   k:= k + 1
   'aux[k]:= (c*'aux[k-1]-'aux[k-i])/(c-1)
   next j
   z:= 'aux[k]
   if abs(z-z0)>1e-8 goto L     " absolute tolerance
   if abs(1-z0/z)>1e-8 goto L   " relative tolerance
call delete('aux)
/E
47
E
#6
:n
elinfk
:z
z = 1/sqrt(1-(k*sin(phi))^2)
/E
E
=r
#1
:r
res1 = Simpson(fun,a,b,n)      " n = 128  yields same precision as  defint1/2
#3
:r
"   If  fun = 'sin                    If  fun = 'elfink
:s
C
#5
:r
"                                              τ b
:s
C
#6
:r
"           τ b                                ³     
:s
C
#7
:r
"           ³                                  ³          d phi
:s
C
#8
:r
"   res1 χ  ³  sin x dx                res2 χ  ³ ΔΔΔΔΔΔΔΔΔΔΔΔΔΔΔΔΔΔΔΔΔΔΔΔ
:s
C
#9
:r
"           ³                                  ³    ___________________
:s
C
#10
:r
"           υ a                                ³   ϋ  1  -  kύ sinύphi
:s
C
#11
:r
"                                              υ a
:s
C
#12
:r
" Integrators with error control:
:s
C
#14
:r
res11 = defint1(fun,a,b)     " Simpson, iterative procedure
#15
:r
res12 = defint2(fun,a,b)     " Simpson, recursive procedure
#16
:r
res13 = defint3(fun,a,b)     " Romberg, iterative procedure
#18
:r
"         ΙΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝ»
:s
C
#19
:r
"         Ί See the explanatory notes in the Function Subsheets Ί
:s
C
#20
:r
"         ΘΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΝΌ
:s
C
=p
=t
=n
=g
#1
:c
.000001
:i
10
:v
Yes
:l
60
:f
66
:m
0
:w
80
:b
No
:a
Yes
:r
Yes
:h
No
:t
1
:n
No
:d
LPT1
%Tr,f,18,1,0,1,0
